function [savings_cv] = generate_cv_results_permutation(sample,pm,split,idx_t,idx_h,rule,reweight,tcost,covariate,winter_prevuse,SMC,savedir,fbase)
% This function estimates EWM rules from the training set and evaluated on the
% test set. 

% Inputs: 
% (1) sample: ross-sectional dataset with electricity consumption, 
% covariates, and propensity scores
% (2) pm: number permutations
% (3) split: share of data in the testing set
% (4) idx_t: index of observations in the training set
% (5) idx_h: index of observations i the testing set
% (6) rule: quadrant, cubic, or univariate rule
% (7) reweight: use density ratio to reweight waves or not
% (8) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (9) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (10) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (11) SMC: =1 use social marginal cost; =0 use retail electricity price
% (12) savedir: directory to save the savings table 
% (13) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Outputs: 
% (1) savings_cv: point estimates of the savings from applying EWM 
% rules derived using the training set on the testing set
%% Initialize training and testing samples
sample_training_set=cell(pm,1);
sample_testing_set=cell(pm,1);

%% verify stratification by opower status and wave is valid using
% cvpartition
for i = 1:pm
sample_training_set{i}=sample(logical(table2array(idx_t(:,i))),:);
sample_testing_set{i}=sample(logical(table2array(idx_h(:,i))),:); 

% verify training set and testing set for each permutation has
% approximately (less than 0.1%) the same share of treated households in RCT 
assert(sum(sample_training_set{i}.opower)/size(sample_training_set{i},1)-...
    sum(sample_testing_set{i}.opower)/size(sample_testing_set{i},1)<0.001)

% verify training set and testing set for each permutation has
% approximately the same share of wave 3, 6, and 7 
    for j=[3,6,7]
        assert(sum(sample_training_set{i}.opower_paper_wave==j)/size(sample_training_set{i},1)-...
            sum(sample_testing_set{i}.opower_paper_wave==j)/size(sample_testing_set{i},1)<0.001)
    end

% verify training set for each permutation has
% approximately the same share of wave 3, 6, and 7 with the pooled sample
    for j=[3,6,7]
        assert(sum(sample_training_set{i}.opower_paper_wave==j)/size(sample_training_set{i},1)-...
            sum(sample.opower_paper_wave==j)/size(sample,1)<0.001)
    end
end

%% run quadrant or cubic rule 
target_wave=0;
sample_wave=0;
savings_cv=cell(pm,1);

if rule=="quadrant"
    perm=1:pm;
    for i=1:pm
        savings_cv{i,1} = generate_results_cross_validation_quadrant(perm(1,i),split,reweight,sample,sample_training_set{i},sample_testing_set{i},sample_wave,target_wave,tcost,covariate,winter_prevuse,SMC,savedir,fbase);
    end
elseif rule=="cubic"
    if pm<100
        perm=pm+1:100;
    else 
        perm=1:pm;
    end    
    
    for i=1:pm
        savings_cv{i,1} = generate_results_cross_validation_cubic(perm(1,i),split,reweight,sample,sample_training_set{i},sample_testing_set{i},sample_wave,target_wave,tcost,covariate,winter_prevuse,SMC,savedir,fbase);
    end   
elseif rule=="univariate"
    perm=1:pm;
    for i=1:pm
        savings_cv{i,1} = generate_results_cross_validation_onedimension(perm(1,i),split,reweight,sample,sample_training_set{i},sample_testing_set{i},sample_wave,target_wave,tcost,covariate,winter_prevuse,SMC,savedir,fbase);
    end       
end
end
